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DUAL  DIFFERENCE  FILTERING:  A  REPLACEMENT 
FOR  INTERPOLATION  AND  SUBTRACTION  TO  DETECT 
CHANGES  IN  MISREGISTERED  SIGNALS 


1.  INTRODUCTION  AND  BACKGROUND 

This  report  introduces  an  alternative  to  interpolation  and  subtraction  as  a  means  of 
detecting  differences  in  a  pair  of  digital  signals.  A  common  application  for  this  method  of 
change  detection  is  to  digitized  imagery  in,  for  example,  renoote  sensing  [1]  or  medical 
imaging  [2,3]  but  it  is  also  useful  in  time  delay  estimation  and  in  many  odier  traditional 
signal  processing  problems  [4].  The  algorithmic  technique  proposed  here  can  relax 
optical/mechanical  alignment  tolerances  for  satellite-based  systems,  reducing  a  vital 
cost/risk  factor  [5].  Improvements  in  the  quality  of  radiological  images  using  these  methods 
can  allow  "reduction  of  injected  contrast  agents  or  X-ray  doses"  [6].  The  new  methods  also 
signiHcantly  enhance  signal-to-clutter  ratios  in  autonomous  surveillance  applications  and  in 
military  search  and  track  systems  [7]  in  which  the  detection  of  moving  targets  is  based  on 
digital  background  subtraction  followed  by  thresholding. 

The  standard  procedure  for  detecting  changes  in  these  applications  is  to  interpolate  the 
first  discrete  signal  at  points  corresponding  to  the  samples  of  the  second,  and  then  subtract 
the  signals.  The  shift  s  between  the  pair  of  sampling  grids  defining  the  discrete  signals  is 
assumed  to  be  known  deterministically  or  from  signal-based  algorithms  [8].  The 
interpolator  applied  to  the  first  signal  usually  takes  the  form  of  a  convolution  with  a  finite- 
extent  kernel  k.  When  evaluated  at  points  shifted  s  from  the  original  grid,  the  discrete  result 
constitutes  a  resampled  digital  signal.  The  second  signal  is  usually  not  resampled. 

In  an  imaging  application,  Stocker  and  Oshagan  [9]  recognized  that  the  spectral 
modification  to  the  first  signal  induced  by  inteipolation  was  lacking  in  the  unaltered  second 
signal.  The  resulting  spectral  mismatch  limited  the  amount  of  clutter  reduction  in  the 
difference  signal.  They  used  a  polynomial  resampler,  the  cubic  B-spline  (CBS),  to  partially 
overcome  this  problem.  CBS  has  tihe  unusual,  noninterpolative  probity  of  not  reproducing 
the  original  discrete  signal  at  the  sample  points,  i.e.,  in  the  limit  s  ->  0.  Therefore, 
resampling  the  second  signal  with  CBS  at  zero  shift  produces  some  spectral  shaping. 
Subtraction  from  the  frrst  signal,  resampled  at  shift  s,  then  achieves  a  marked  reduction  in 
residual  clutter. 

One  way  to  eliminate  spectral  mismatch  completely  is  to  apply  a  standard  symmetric 
interpolator  to  the  first  signal  at  the  shift  s/2,  and  to  the  second  at  -  s/2.  However,  better 
performance  can  be  realiz^  if  the  restrictive  idea  of  interpolating  is  discarded.  The  goal  of 
interpolation  is  to  reproduce  an  underlying  continuous  signal.  It  is,  therefore,  appropriate 
when  a  comparison  is  made  to  an  unaltered  second  signal.  But  if  the  actual  goal  is  change 
detection,  this  approach  is  unnecessarily  limiting. 

This  report  considers  the  general  effects  of  filtering  both  signals  prior  to  subtraction 
and  develops  a  mathematical  method  for  evaluating  and  optimizing  these  Dual  Difference 
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Filters  (DDFs).  Ccmventional  interpolation  is  seen  as  a  special  case,  when  the  sectmd  signal 
is  n(A  filtered  Interpolators  with  nontrivial  behavior  at  zero  shift  like  the  cubic  B-spline  are 
seen  as  suboptimal  types  of  DDF. 

2.  ANALYSIS 

With/fxj  representing  the  underlying  signal  intensity  at  position  x,  sample  values  of 
(where  n  is  an  integer)  constitute  the  digital  signal,  which  may  convolved 

with  a  kernel  to  produce  a  discrete  sequence 


(1) 

n  * 

In  practical  applications,  this  sum  is  finite  because  is  nonzero  for  only  a  few  values 
of  m  near  m  =  0.  The  dependence  of  k\  on  s  arises  because  the  choice  of  filter  generally 
depends  on  the  separation  of  the  sampling  grids. 

Similarly,  the  second  digital  signal /(n  +  s),  which  is  a  sampled  version  of /  at  the 
shifted  points,  can  be  filtered  to  produce  a  second  sequence,  using  a  (possibly)  different 

kernel 


l4**«n-«)/(n  +  s).  (2) 

n 

It  is  convenient  to  redefine  ki  to  allow  for  continuous  arguments 

(3) 

where  n  is  an  integer.  When  Eq.  (3)  is  substituted  into  expression  (1),  the  filtered  sequence 
in  (1)  assumes  the  form  of  a  convolution  evaluated  at  the  point  m  +  s,  which  is  already  the 
form  in  (2).  The  sequences  in  expressions  (1)  and  (2)  now  can  be  treated  as  sample  vdues 
of  continuous  functions /i,/2,  where 


fi(x)=I,ki(x-n)f(n) 

n 


(4a) 


and 

-s)f{n  +  s) ,  (4b) 

n 

that  are  evaluated  at  x  =  m  +  5.  Equation  (4a)  is  a  standard  form  for  an  interpolation 
function /i  produced  by  a  kernel  ki. 

We  will  usually  assume  that  0  ^  ^  1/2.  That  is,  the  second  signal  is  shifted  to  the 
right  of  the  first  by  no  more  than  half  a  sample.  This  can  always  be  arranged  by  relabelling 
the  integer  sample  points  and/or  interchanging  the  order  of  the  signals,  if  necessary.  Also, 
the  treatment  here  is  in  one  dimension  only.  Higher-dimensional  filtering  can  be 
accomplished  by  a  cascade  of  two  one-dimensional  operations  or,  more  generally,  by 
treating  the  arguments  x,  n,  and  s  as  vectors. 
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In  conventional  change  detection,  before  a  subtraction  is  made,  one  of  the  signals — 
say,  the  first — is  interpola^,  and  tiie  other  is  unaltered.  For  this  meAod,  a  good  kernel  k\ 
produces  an/i  that  is  close  to/,  and  the  second  kernel  is  trivial: 

5w.  (5) 

where  5  is  the  Kronecker  delta.  However,  the  requirement  that  be  a  good  interpolator  is 
too  restrictive  if  the  real  goal  is  change  detection.  Nevertheless,  a  review  of  the 
interpolation  application  is  instructive. 

Because  the  sampling  grids  {«}  and  {n  +  5}  are  defined  here  as  having  unit  spacing, 
the  Nyquist  frequency  is  VNyq  =  1/2  cycles/sample.  If  the  Fourier  transform  of/ has  no 
components  with  frequency  v  ^  1/2  ,/is  called  "oversampled."  According  to  the  Nyquist 
Reconstruction  Theorem,  the  kernel 


ki{x)  =  sinc(j:)  =  (6) 

may  be  used  in  Eq.  (4a)  to  reproduce  /exactly  (i.e.,/i  =/),  whenever/ is  oversampled. 

In  practice,  k\  must  be  of  finite  support  (the  range  of  x  over  which  k\  is  nonzero)  to 
make  expression  (1)  or  Eq.  (4a)  contain  only  a  finite  number  of  terms.  For  example,  the 
sine  function  can  be  truncated  at  ±  iV/2  to  make  an  N-point  interpolator.  Two  other 
common  interpolators  are  nearest  neighbor  (NN)  and  linear  (LIN).  Figure  1  illustrates 
k\{x)  for  NN,  LIN,  and  SINC.  The  supports  for  NN  and  LIN  are  1  and  2,  respectively. 
The  SINC  has  infinite  support,  but  when  tmneated  at  ±  N/2,  its  support  or  order  becomes 
N. 
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The  method  of  analysis  developed  here  relies  heavily  on  Fourier  transform  theory. 

The  transform  operator  is  denoted  by  T,  and  upper  and  lower  cases  denote  transform  pairs. 
Also,  the  symbol  indicates  a  correspondence  between  real  and  transform  space:  for 
example, 

(ffc)  (v)  s  A:(v)  s  j  e-ifox  k(x)  dx  (©  s  27cv)  (7) 

means  the  same  thing  as 

k^K.  (8) 

The  shift  and  convolution  theorems  are  then 

g(pc  +  s)  <z>e2J«sv  C(v) 

e-2i«xvo  g(x)  <=>  G(v  +  Vo) ,  (9) 

and 

g<x)A(r)«(G*H)(v) 

{g*h)Oc)^CHy)H(v).  (10) 

Another  inqxntant  thecnem  [10,1 1]  relates  the  comb  function,  an  infinite  sum  of  Dirac 


delta  functions. 

comb(x:)s  j;5(r-n) , 

(11) 

to  its  transform 

comb(x)  <=>  CX3MB(v) . 

(12a) 

Equations  (1 1),  (12a),  and  (7)  immediately  imply  that 

COMB(v)  =  Zc-2’w"v 

(12b) 

n 


Another  important  relation  that  follows  directly  from  the  definitions  of  the  continuous 
convolution  operation  *  and  the  comb  is 

XgOc -  n) h(n)  =[g  ♦  (comb  •  /i)]Cc) ,  (13) 

n 

which  is  valid  for  arbitrary  g,  h.  A  special  case  of  Eq.  (13)  is  also  useful: 

2  gCt  -  n)  =  (^  *  comb)(x) .  (14) 

n 

These  relations  can  be  used  to  express  a  discrete  sum  of  squared  sample  values  of  a 
function  g  in  terms  of  its  Fourier  transform  C.  Letting  g  g^,x—*0,  and  n  — >  -n  in  Eq. 
(14)  results  in 


I«2(n)=(g2  ♦  comb){0) . 


(15) 
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Next,  successively  using  Eqs.  (10)  and  (12a),  the  definition  of  and  Eq.  (10)  again  in 
Eq.  (15)  produces 

.comb))](0)  =  [r‘(r(*2)COMB)](0) 

=  j  <<v'[r(g2)(v-)COMB(v')]=  f(iv'[G.C](v')CX3MB(v').  (16) 

The  frequency-space  version  of  Eq.  (1 1)  is  now  used  in  Eq.  (16)  along  with  the  definition 
of  convolution: 


=  dv'5(v'-n)| dvG(v'-v)  G{v) 


(17) 


dvG{/t-v)  G{v). 

We  assume  for  convenience  that  g  is  real.  Then 

G(v)  =  G*(-v)  Vv, 


(18) 


and  from  Eq.  (17), 


=  I  dvG*(v-n)  G(v). 


(19) 


If,  furthermore,  the  sequence  g(n)  results  from  oversampling  the  continuous  function 
g,  so  that  G(v)  =  0  for  I V  I  ^  1/2 ,  then  Eq.  (19)  simplifies  to 


Ig2(/»)=  dv|G(v)|2.  (20) 


which  is  closely  related  to  Parseval's  Theorem,  with  the  usual  roles  of  position  and 
frequency  interchanged.  Note  the  distinction,  however:  in  Eq.  (20),  G  is  the  Fourier 
transform  of  the  underlying  function  g,  not  of  the  Fourier  series  {g(n)}.  These  two 
transforms  are  equal  if  g  is  oversampled,  but  even  if  it  is  not,  Eq.  (20),  which  we  will  call 
the  discrete  Parseval's  'Aeorem,  still  holds  in  an  average  sense  describe  below. 


Note  from  Eq.  (20)  that  for  any  oversampled  g,  Z  gHn)  is  indepnendcnt  of  where  the 
sampling  grid  is  laid  down  relative  to  the  continuous  signal  g(x).  Shifting  the  grid  by  an 
amount  t  in  one  direction  is  equivalent  to  shifting  g(x)  in  the  otiier,  i.c.. 


g(j:)^  g(jc  +  /) 

which,  according  to  Eq.  (9),  in  Fourier  space  is  equivalent  to 

G(v)-»  G{v)e2Jtity', 


(21) 

(22) 
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under  which  change  the  right  side  of  Eq.  (20)  is  invariant.  Because  the  value  of  Z  gHn)  is 
indej^ndent  of  grid  placement  for  an  oversampled  g,  we  call  Eq.  (20)  the  strong  form  of 
the  discrete  Parseval's  Theorem.  When  g  is  not  oversampled,  a  weaker  version  still  holds, 
in  which  g  is  randomly  shifted  with  respect  to  the  sampling  grid. 

In  particular,  under  the  translation  of  the  function  g  by  the  amount  t  e  [0,1],  according 
to  Eq.  (22),  Eq.  (19)  becomes 


^  glititn 
n 


I 


dvG*(v-/i)  G(v). 


(23) 


Averaging  Eq.  (23)  over  /;  0  ->  1  yields  Eq.  (20).  That  is,  the  weak  form  of  the  discrete 
Parseval's  Theorem  is  also  express^  by  Eq.  (20),  if  the  left-hand  side  is  understood  as  an 
average  over  all  placements  of  the  sampling  grid. 


Methods  similar  to  those  above  can  also  be  used  to  generalize  both  forms  of  the 
discrete  Parseval's  Theorem  and  to  derive  a  simple  formula  characterizing  the  performance 
of  any  DDF.  To  compute  the  differencing  error  associated  with  the  DDF  (Jfci,  k^),  we  let 


g(»)  (x)  =/i  (x  +  s)  (x  +  s)  (24) 

(see  Eq.  (4)).  Then  the  squared  error  in  the  difference  signal  is 

d}='L\fi(n+s)-f<^Kn  +s)P,  (25) 

n 

or 

d?  =  I[sW(i.lP.  (26) 

Next,  finom  Eqs.  (24)  and  (9), 

G(^)  (v)  =  [  Fi  (v)  -  (v)]  (27) 

The  asymmetry  in  the  (5) -dependence  introduced  in  Eqs.  (3)  and  (4)  accounts  for  the 
corresponding  asymmetry  in  the  F,-  of  Eq.  (27).  This  notation  was  introduced  to  facilitate 
use  of  the  discrete  Parseval's  Theorem. 


If  the  sampling  grid  were  dense  enough  that  g^^Xn)  represented  an  oversampling  of 
g^^Xx),  then  Eq.  (20)  could  be  applied  immediately  to  Eqs.  (26)  and  (27),  with  g  -> 

But  the  functions /•  deOning  gW  in  Eq.  (24)  arc  usually  not  bandlimited,  and  so  neither  is 
because  the  transformations  of  Eq.  (4)  introduce  aliasing  whenever  either  kernel  ki  is 
of  finite  support.  Therefore,  g  cannot  be  oversamplcd,  even  if  the are,  and  the  strong 
Parseval’s  Theorem  does  not  apply. 

Nevertheless,  when  Eqs.  (24),  (4a),  and  (4b)  define  g,  Eq.  (20)  still  holds  whenever 
/is  oversampled,  even  if  g  is  not.  This  key  property  will  permit  the  simple  characterization 
of  the  error  (Eq.  (25))  as  a  power  spectrum  that  is  the  product  of  an  attenuation  factor  and 
the  power  spectrum  of  /.  The  attenuation  factor  depends  only  on  the  kernels  ki  and  the 
shift  s. 
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Using  Eq.  (13)  in  Eq.  (4a)  with  h  =  /and  g  =  k\,  taking  the  transform,  and  applying 
Eq.  (10)  results  in 


Fi  (v)  =  Kx  (v) [CXDMB  ♦  F]  (v) ,  (28) 

which  is  usually  not  bandlimited  because  neither  factor  is,  as  will  be  demonstrated. 
Consequently,  G  in  Eq.  (27)  is  usually  not  bandlimited. 

The  frequency  version  of  Eq.  (14)  allows  us  to  rewrite  Eq.  (28)  as 

Fi(v)  =  Ari(v)IF(v-m).  (29) 

m 


By  similar  means,  it  can  be  shown  that 

(v)  =  (v)  (30) 

m 

Even  if  F  is  bandlimited,  the  second  factor  in  Eq.  (29)  or  Eq.  (30)  is  nonzero  for  arbitrarily 
large  v.  So  is  the  first  whenever  the  kernel  is  of  finite  support,  i.e.,  in  all  practical 
cases.  Consequently,  g  is  not  bandlimited,  as  claimed,  because  the/  are  not. 

Nevertheless,  using  Eq.  (26)  along  with  Eqs.  (19),  (27),  (29),  and  (30)  results  in 

4  =  X  I  d\  XF*(v-n-m)[^i(v-n)-e'2’i«"‘^:^^^(v-n)]*  x 

n  \  m 

e2nisn^f{y.j)  [A:i(v)-e*2««;4*nv)] ,  (31) 

j 

which  we  now  show  to  be  independent  of  grid  placement  whenever/ is  oversampled. 
Assuming  then  that 

F(v)  =  0  for  |vl  ^  (32) 

we  examine  a  typical  cross-term  of  Eq.  (31): 

_  £  e27t««  I  dv  [^1  (v  -  n)]*  X  fV  -  n  -  m)  4*\v)  X  ^  (v  -  j)  .  (33) 

n  J  m  j 

Because  of  Eq.  (32),  the  integration  in  1^.  (33)  kills  all  summand  terms  except  when  n  +  m 
=  j.  Notice  also  that  averaging  over  grid  placement  t,  as  described  in  Eqs.  (21)  and  (22), 
has  the  same  effect.  This  means  that  the  following  results  hold  in  the  weak  sense  also,  for 
which  Eq.  (32)  need  not  apply. 
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Equation  (33)  thus  reduces  to 

dv [Ki  (v  - n)]*  S I F(v - n - m) P (v)  .  (34) 

m 

Changing  variables,  v-»v  +  n  +  /nin  Eq.  (34)  produces 

- 1  dv  lF(v)  1 2  X  (v  +  m)]*  4*^  (V  +  n  +  m)| .  (35) 

J  njn  I 


es, 

-/ 


The  other  cross-terms  of  Eq.  (31)  produce  results  similar  to  expression  (35),  but  with 
different  terms  within  the  braces.  The  final  result  may  be  written  as 


=||F{v)|2|£, 


(v)l2dv  , 


(36) 


which  includes  an  attenuation  factor  that  is  the  squared  magnitude  of  a  "complex  error 
factor"  £,(v): 


Es  (v)  s  Ki  (v  +  «)  -  4'^  (v  +  n)] .  (37) 

n 

The  quantity  IF(v)|2|£,(v)|2  will  be  called  the  error  spectrum.  The  error  factor  ^^(v) 
depends  only  on  the  kernels,  and  not  on  the  signals  being  processed.  The  independence  of 

Eq.  (36)  from  the  phase  of  F  means  that  dj  does  not  depend  on  the  placement  of  the 

sampling  grids  when  the  signal  is  oversampled,  even  though  the  individual  terms  of  Eq. 
(25)  do.  Furthermore,  as  indicated  earlier,  Eq.  (36)  holds  also  for  undersampled  signals 
when  averaged  over  dl  such  placements.  Therefore,  the  result  expressed  by  Eq.  (36)  has 
strong  and  weak  forms,  just  as  the  discrete  Parseval's  Theorem  has  (Eq.  (20)). 

The  formula  in  Eq.  (37)  is  most  useful  when  the  sum  contains  only  a  finite  number  of 
terms,  that  is  when  the  Fourier  transforms  of  the  kernels  are  of  finite  support  For  example, 
for  perfect  (inband)  interpolation  ii(jc)  =  sinc(x),  and  ^ri(v)  has  the  value  1  inband  (-1/2  S 

V  ^  1/2),  and  is  zero  out  of  band.  If  the  second  signal  is  left  unaltered,  then  k2\x)  can  be 

chosen  to  be  any  function  satisfying  Eq.  (5);  a  convenient  choice  is  =  sinc(jr).  Then 
Ar2(v)  =  Ari(v),  and  Eq.  (37)  can  be  used  to  show  that  for  sine  interpolation  the  error  factor 
is  given  by 


[£5  (v)  1  =  2 1  sin  |7cs  int |v  *0 1 .  (38) 

where  "int"  means  "largest  integer  less  than  or  equal  to."  (For  example,  int[-.5]  =  -1.) 
Thus,  the  error  factor  is  zero  inband,  but  not  out  of  band,  where  many  other  interpolators 
are  superior  [12]. 

In  many  practical  situations,  it  is  the  real-space  form  of  a  kernel  that  is  of  finite 
support,  and  so  the  next  goal  is  to  express  Eq.  (37)  in  terms  of  ki  instead  of  ATj. 
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The  first  term  of  Eq.  (37)  is,  by  definition, 


iki{x)e-2wXv+nWjc. 


(39) 


The  real- space  version  of  Eq.  (12b)  along  with  Eq.  (1 1)  then  permits  us  to  rewrite  Eq.  (39) 
as 


j  ki  (x)  comb  (s  -  x)  dx  =  X  e-2niv(5-w)  ki  (s  -i-  /i) .  (40) 

A  similar  expression  for  the  second  term  of  Eq.  (37)  leads  to 

£s{v)  =  X  ki  (5-i-/i)-4^^(n)]  =  X  (/i) - ,  (41) 

n  n 

where  Eq.  (3)  has  also  been  used.  Equation  (41)  is  an  alternate  form  for  Eq.  (37);  it 
expresses  the  error  factor  directly  in  terms  of  the  kernels  instead  of  their  Fourier 
transforms. 


According  to  expressions  (1)  and  (2),  the  values  k^i\n)  and  l^\n)  are  tap  weights  of 
digital  filters  operating  on  the  pair  of  discrete  signals  f{n)  and/(n  +  s).  Calling  the  transfer 
functions  of  these  filters  and  \  i.e.. 


(42a) 

(42b) 

the  complex  error  factor  may  be  written  compactly  as 

E,  (v)  =  //f  ^(v)-  (v) 

(43) 

Note  that  the  sums  in  Eqs.  (42a)  and  (42b)  are  finite  in  practical  applications.  In  the 
conventional  method  of  change  detection,  in  which  the  first  signal  is  interpolated  and  the 
second  is  unaltered,  Eq.  (5)  applies,  and  the  general  error  factor  for 
interpolation/subtraction  reduces  to 

£,(v)  =  e-»w^//</>(v)- 1  . 

(44) 

The  form  for  interpolation  error  that  results  from  inserting  Eq.  (44)  into  Eq.  (36)  is 
equivalent  to  one  derived  in  Ref.  13  for  the  special  case  of  a  purely  sinusoidal  input  signal 
with  a  sub-Nyquist  frequency. 
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As  an  (  xample,  it  can  be  shown  that  the  perfect  inband  interpolator  (Eq.  (6))  has 
transfer  funcion 


(sine  interpolation) .  (45) 

Equation  (38),  which  was  derived  from  the  continuous-transform  representation  in  Eq. 
(37),  can  be  recovered  by  substituting  Eq.  (45)  into  Eq.  (44). 

Equations  (36),  (37),  and  (41)  or  (43)  express  the  resampling  error  at  any  fixed  shift 
s  between  a  pair  of  sampling  grids  in  terms  of  the  spectrum  of  the  underlying  continuous 
signal  and  an  error  factor  that  depends  only  on  the  method  of  filtering  in  expressions  (1) 
and  (2).  Furthermore,  if  /  is  not  oversampled,  then  Eq.  (36),  which  constitutes  a 
fundamental  theorem  of  resampling  error,  still  holds  in  the  weak  sense,  i.e.,  averaged  over 
all  possible  placements  of  the  two  sampling  grids  relative  to  the  signal.  Finally,  if  the 
function /is  reg^ded  as  one  realization  of  a  stationary  stochastic  process,  then  taking  the 
ensemble  average  of  Eq.  (36)  yields  a  third  form  of  the  theorem,  which  relates  the  mean 
squared  error  to  \E\^  and  P(\)  =  <IF’(v)l^>e^n,bie’  which  is  the  power  spectrum  of  the 
stochastic  process. 

3.  OPTIMAL  DUAL  DIFFERENCE  FILTERS 

The  three  forms  of  the  fundamental  theorem  can  be  used  to  design  DDEs  satisfying  a 
variety  of  criteria,  such  as  error  minimization.  For  example,  Eqs.  (36)  and  (41)  may  be 
combined  to  give 

-2  £  R{n  +s)ki{n+m+s) f^^(ni)  (46) 

ntjn 

with 


/?(x)s  cos(x©)|F(v)pdv.  (47) 


For  the  interpolation  application  (see  Eq.  (5)),  this  reduces  to 

d}=  R{Qi)-2'ZR(n  +  s)k\(n  +  s)+  F{n)k\{n  +  m  +  s)ki{m  +  s) .  (48) 

n  m/i 

Minimization  of  in  Eq.  (48)  with  respect  to  the  variables  k\{n  +  s)  results  in  an 
equation  for  the  optimal  N-point  interpolation  kernel 

£/?  (/n)  («  +  5)  =  /? (n  +  j) ,  (49) 

m 

valid  for  those  n  for  which  kiin  +  s)  in  Eq.  (48)  is  nonzero,  i.e.,  for  the  N  values  of 
n:  int[-(N-  l)/2] ,...,  int[(A-  l)/2].  Comparing  Eq.  (49)  with  Eq.  (4a)  shows  that  the 
former  is  equivalent  to  the  requirement  that  ki  serve  as  a  perfect  interpolator  for  the 
function  R(x)  for  x  e  [-N/2,  A/21,  i  e-.  over  the  support  defining  the  order  of  ki.  Equation 
(49)  has  appeared  before  [1],  but  only  as  the  condition  minimizing  an  ensemble 
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mean- squared  error  for  a  stationary  stochastic  process,  that  is,  in  a  more  limited  application 
than  shown  here. 

Because  Rim)  is  symmetric,  Eq.  (49)  is  a  Toplitz  equation,  for  which  general  iterative 
solution  procedures  have  been  developed  [14,15].  Thus,  in  principle,  kiis  +  n)  can  always 
be  found  as  a  function  of  2N  numbers,  R(m)  and  Rin  +  s),  with  m  =  0,...,  -  1  and 

n  =  int[-(A/  -  l)/2] ,...,  int[<(iV  -  l)/2].  Solutions  to  (49)  were  found  in  Ref.  12  for 
various  power  spectral  models:  uniform,  Gaussian,  power  law,  and  Lorentzian. 

Error  minimization  equations  can  be  derived  for  the  mote  general  case  of  DDEs  in  Eq. 
(46),  but  now  at  least  one  extra  constraint  must  be  added,  lest  the  trivial  solution  k\in  +  s) 

=  k2\n)  =  0  result.  The  nature  of  such  constraints  depends  on  that  of  the  change  that  is  to 

be  detected.  It  may  be  an  extensive  change,  such  as  in  a  region  of  blood  flow  on  an 
angiogram,  or  large-scale  agricultural  evolution  as  seen  from  Earth  orbit.  It  may  be  a  local 
change,  such  as  tire  motion  of  a  small  object.  The  general  goal  is  to  combine  the  two 
signds  after  Altering  so  that  only  such  changes  persist. 

In  dual  difference  filtering,  at  least  one  of  the  individual  filters  ki  must  be  designed  to 
preserve  the  change  that  is  to  be  detected.  As  an  example,  for  extensive  changes,  an 
appropriate  choice  of  constraint  is  to  require  perfect  dc  response  of  each  filter.  This  will 
tend  to  preserve  local  averages.  Similarly,  for  detection  of  small  moving  targets  of  known 
shape,  the  filters  may  be  constrained  to  conserve  target  energy  or  peak  amplitude. 

Generally,  the  inclusion  of  constraints  in  the  DDF  optimization  problem  can  be 
achieved  by  intnxlucing  Lagrangian  multipliers  Xi  into  Eq.  (46): 

Z  R{n)[ki{n  +  m+  s) ki{m  +  s)  +  +  m) (m^ 

mjn 

-2  Rin  +  s)ki{n+m  +  5)  (m) - 2 S Xi +  s),  4*\n)))  •  (50) 

m/i  i 

The  functions  c,-  represent  constraints  on  the  digital  filters.  The  error  minimization 
equations  then  result  from  differentiating  Eq.  (50)  with  respect  to  the  tap  weights  kiin  +  s) 

and  k^\n),  and  X,,- : 


X/?  M  Jtiln  -m  +  s)-'ZR{m  +  s)  iSpin -m)  =  £X,,- 

mm  I 


dcj 

dki{n+s) 


(51a) 


and 


X  R  {m)  l^\n  -m\-JiR{m+s)  ki{n  +  m  +  s)  =2  X-i  — — , 
m  m  I  dk2  \n) 


(51b) 

(51c) 


In  many  optimization  problems,  the  solutions  for  ku  fp  show  a  symmetry  reflecting 
the  fact  that  neither  of  the  two  original  discrete  signals  is  preferred.  The  symmetry  is 


k\{n  +  s)  =  4^^  (  -m) , 


(52a) 
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or,  from  Eq.  (3), 

k^f^(n)  =  I^H-n),  (52b) 

which  is  illustrated  in  Fig.  2  for  two-  and  three-point  filters.  In  the  remainder  of  this  report, 
the  symmetry  in  Eq.  (52)  will  be  assumed  for  all  DDFs,  and  the  subscript  on  ki  will  be 
omitted.  Other  symmetries  are  considered  in  Ref.  16. 


r(n) 


I  I 


k(s)  k(s-l) 


f(n) 


f(n-i ) 


FRAME  I 


- FRAME  2 


f(n*l) 


f  A  ^  ^  f  f 


Three-point  DDF 


k(s*1)  k(s)  k(s-l) 

Fig.  2  —  The  usual  symmetry  in  DDF  filters.  The 
double-arrow  lines  connect  those  sample  points  in  the 
two  digital  signals  that  are  multiplied  by  equal-value 
tap  weights  from  the  corresponding  kernels. 
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For  the  symmetric  case,  the  DDF  error  factor  (Eq.  (41))  reduces  to 

£,(v)  2  X*{-s  +  «)sin[^©|/i  +  ,  (53) 

in  which  an  irrelevant  unimodular  phase  factor  has  been  omitted.  Also,  the  residual  error 
after  subtraction  (Eq.  (46))  simplifies  to 

4  =  2  1  {R{n)k{n  +  m  +  s)k{m  +  s)  —R{n  +  s)k{n  +  m  +  s)k{—m  +  s)} ,  (54) 

m,n 

and  now  there  are  fewer  independent  optimization  equations: 

'ZR  {m)  k  {n  -  m  +  s)  -  ZR  im  +  s)k{m-n  +  s)  =  Z^i  \  (55a) 

m  m  i  o 

Ci{[k{n  +  s)))  =  0  (55b) 

than  in  the  general  case  (Eq.  (51)). 

The  only  type  of  constraint  considered  here  in  detail  is  that  of  a  single  linear  one, 
which  can  be  written  as 


c  ({/:  (rt  +  5)})  s  5^  (n  +  s)  - 1  =  0  .  (56) 

n 

If  r_„=  1,  Eq.  (56)  describes  the  "dc"  constraint.  Then  the  transfer  functions  in  Eq.  (42) 
have  unit  response  at  v  =  0.  Alternatively,  (fn)  may  describe  the  amplitude  of  a  target 
image  whose  interframe  motion  is  to  be  detect^  by  use  of  a  DDF.  Gaussian  targets  will  be 
considered,  with  a  standard  deviation  of  a  samples  and  unit  peak  value 

r„  =  e-n2/2a2  (57) 

In  all  cases  considered,  to  has  been  normalized  to  unity.  The  values  o  =  1/2, 1  will  be  used 
as  examples.  In  imaging  applications  with  matched  optics  design,  these  correspond 
respectively  to  undersampled  and  nearly  Nyquist-sampled  images  of  point  objects  [12]. 

The  constraint  in  Eq.  (56)  then  means  that  the  peak  signal  value  is  reproduced  by  the 
filter  (if  the  peak  target  value  lies  on  the  first  sampling  grid).  It  often  happens  that  an 
asymmetrical  form  of  the  tap  weights  [k^in)]  transforms  a  symmetrical  target  into  an 
asymmetrical  one,  with  a  consequent  shift  in  the  peak  location.  The  filtered  pe^  can  then 
be  larger  than  the  unfiltered  one.  In  any  case,  Eq.  (56)  requires  that  the  peak  value  in  the 
filtered  signal  be  no  less  than  in  the  unfiltered. 

In  all  DDF  examples  studied  here,  error  minimization  requirements  are  satisfied 
simultaneously  with  a  constraint  that  peak  signal  be  conserved.  These  are  always  equivalent 
to  a  signal-to-clutter  criterion.  For  example,  requiring  that  the  peak  signal  be  doubled 
instead  of  conserved  generally  also  doubles  the  filter  tap  weights  and,  from  Eq.  (46)  or 

(54),  the  rms  clutter  as  well. 
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With  the  constraint  in  Eq.  (56),  Eq.  (55a)  becomes 

'^R{ni)k(n-m  +  s)-'^R{m+s)k(m-n+s)  =  Xt.n ,  (58) 

m  m 

valid  for  n  =  int[-(iV  -  l)/2] int[(N  -  l)/2].  These  linear  equations  in  N  variables 
k(n  +  s)  should  be  compared  to  Eq.  (49),  the  corresponding  optimization  equations  for  the 
interpolation  problem.  The  explicit  constraints  in  Eq.  (56)  are  replaced  there  by  a  more 
severe  one — ^the  second  signal  is  left  unfiltoed. 

Before  applying  these  results  to  specific  problems,  we  note  that  as  long  as  the 
quant’^ies  R(m)  and  R(m  +  s)  are  known,  Eq.  (58)  can  be  used  to  find  the  optimal  kin  +  s) 
even  if  the  value  of  s  is  unknown,  because  all  the  s-dependence  in  Eq.  (58)  is  implicit, 
appearing  in  the  arguments  of  functions.  This  means  that  change  detection  might  be 
possible  without  knowledge  of  the  shift  s  between  sampling  grids,  if  estimates  of  Rim)  and 
Rim  +  s)  are  available.  And  these  may  be  provided  by  interpreting  /?  as  an  autocorrelation 
function  and  then  constructing  periodograms  [17];  (a)  from  either  signal,  to  estimate  Rim), 
and  (b)  from  both  to  estimate  Rim  +  s).  This  technique  introduces  errors  by  using 
estimates  of  R,  but  it  completely  eliminates  errors  from  a  prior  step  that  we  have  ignored: 
estimation  of  the  shift  s.  The  competition  between  these  errors  depends  on  the  performance 
of  registration  methods  and  is  not  studied  here. 

3.1  Power>Law  Spectrum 

A  common  analytic  model  for  image  power  spectra  is  a  low  inverse  power  p  of 
frequency,  typically  with  p  »>  2.  We  will  examine  optir^  DDEs  f«  the  spectrum 

|F(v)|2-1/v2.  (59) 

Besides  its  relevance  to  imagery,  this  choice  is  motivated  by  the  simplicity  of  the  optimal 
interpolator  solution,  which  is  known  to  all  orders  and  can  therefore  be  compared  easily 
with  the  DDF  solutions.  Also,  the  method  used  here  to  handle  the  dc  singularity  is 
generalizable  to  other  singular  spectra.  In  practice,  at  some  cutoff  the  low-frequency  form 
of  Eq.  (59)  becomes  bounded,  keeping  the  real-space  image  variance  finite.  The  following 
analysis  shows  that  the  relative  performance  of  most  resamplers  is  insensitive  to  the  cutoff. 

We  will  treat  Eq.  (59)  as  the  limiting  case  of  the  spectrum 

|P{v)P-^J-^  (60) 

V2  +e2 

as  e  ->  0.  The  function  R  in  Eq.  (47)  is  then  easily  evaluated: 

/?(x)  f  pW 

where 

pse-2ne  (e^O).  (61) 
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For  two-point  filtering,  Eq.  (58)  reduces  to 


1 

whose  solution  is 


R{\)-R{l-s)  -1 

R(0)-R{2-s)  -ti 

h  0 

k  (i)  =  m3i  /  det 
k(s-\)=-myil  det 


k(s) 

o' 

k  (5-1) 

= 

0 

X 

1 

(62) 


(63) 


with  niij  the  cofactors  of  the  matrix  in  Eq.  (62)  and  det  the  determinant.  In  particular, 

m3i  =R{0)-R{2-s)-ti[R(l)-R{l-s)] 
m22  =^R{l)-R{l-s)-ti[R(0)-R(s)] 
det  =m3i  -/i  m32  .  (64) 


The  limit  e  ->  0  corresponds  to  p  -»  1,  in  which  case  (e/n)/f(jc)  1  V  x;  Eq.  (64) 

then  shows  that  Eq.  (63)  is  indeterminate.  Application  of  IHdpital's  Rule  yields 

k{s)  =  [s{\-tx)-2\l  det 

^(jf-l)  =  -5(l  +ri)/tfer 

tfef  =  2(s-l)-5(i  +  r,f,  (65) 

which  describes  the  optimal  two-point  DDF  for  a  p  =  2  power-law  spectrum.  Except  when 
t\  assumes  one  of  the  two  values  -1  ±  V2,  this  kernel  is  a  rational  function  of  s,  rather  than 
a  polynomial.  Note  that  the  optimal  interpolation  kernels  for  power-law  spectra  are 
polynomials  [12],  as  are  most  interpolators  in  popular  use.  Thus,  the  DDF  approach 
inti^uces  a  new  class  of  kernels  as  a  tool  for  change  detection. 

For  the  power  spectrum  of  Eq.  (59),  the  minimum-error  interpolator,  i.e.,  the 
solution  to  Eq.  (49),  is  LIN  for  all  N  2:  2  [1].  TTiis  is  a  peculiarity  of  the  p  =  2  spectrum, 
that  the  optimal  interpolator  has  support  2  instead  of  N.  Figure  3  plots  the  error  factors 
IEj(v)l  for  each  of  three  interpolators;  the  optimal  (LIN),  Cubic  Convolution  [18],  which  is 
a  popular  four-point  interpolator,  and  the  Low-Frequency  optimal  four-point  interpolator 
(LF-4),  which  is  designed  to  have  the  smallest  possible  error  factor  in  the  vicinity  of  v  =  0 
[12]  and  is,  therefore,  appropriate  for  the  singular  spectrum  in  Eq.  (59).  At  s  =  0,  the 
errors  are  zero;  at  s  =  1/2,  they  are  usually  the  worst.  Our  examples  are  always  at  the 
intermediate  shift  s  =  1/4,  which  choice  explains  (see  Eq.  (41))  the  Av  =  4  periodicities  in 
the  error  factors. 
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FREQUENCY  (CYCLES/SAMPLE) 

Fig.  3  —  Emv  factors  Ifsfv)!  at  the  shift  s  s  0.2S  for  three 
interpolators:  Linear,  Cubic  Convolution,  and  LF-4;  and  one 
two-point  DDF:  the  optimal  for  a  p  =  2  power  law  with  the 
dc  constraint  The  Nyquist  frequency  is  0.5  cycles/sample 


Figure  3  also  shows  the  error  factor  for  the  optimal  two-point  DDFs  in  Eq.  (65), 
using  the  dc  constraint  (t„  =  1  in  Eq.  (56)).  Rgure  4  shows  the  corresponding  error  spectra 
(the  integrand  of  Eq.  (36))  for  the  case  of  a  p  =  2  signal  spectrum,  and  it  lists  the  relative 
areas  under  these  curves.  By  Eq.  (36),  these  are  also  the  relative  mean  squared  errors  in  the 
filtered  difference  signals.  Note  that  the  optimal  DDF,  even  for  N  =  2,  ouq)erforms  all 
interpolators  regardless  of  their  order,  because  it  outperforms  the  optimal,  LIN.  The 
polynomial  interpolators  (Cubic  Convolution  and  LF-4)  are  included  to  illustrate  their 
superior  performance  below  the  Nyquist  frequency  (v  =  1/2),  an  important  point  in  our 
later  discussions. 


FPEOUENCi  (CYCLES/SAMPLE) 


Fig.  4  —  Error  spectra  l£g(v)  F(v)|2  when  the  image  spectrum 
is  a  p  «  2  power  law,  for  three  interpolators  and  one  DDF  at 
the  shift  s  «  0.2S 
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Figures  5  and  6  show  the  effects  of  imposing  on  individual  DDF  filter  performance 
the  three  types  of  constraint  discussed  earlier.  Figures  7  and  8  show  the  corresponding 
results  for  the  p  =  2  optimal  kernels  when  N  =  4.  As  always,  a  shift  s  =  1/4  is  assumed. 
The  increasing  errors  with  decreasing  o  occur  because  smaller  values  of  the  standard 
deviation  correspond  to  targets  with  a  larger  high-frequency  content  Therefore,  conserving 
the  peak  signal  value  competes  increasingly  with  minimizing  difference  clutter,  for  which 
low  frequencies  are  important  for  power-law  spectra.  The  dc  constraint  may  be  considered 
a  limiting  version  of  the  Gaussians  as  o  ->  «>. 


FREQUENCY  (CYCLES/SAMPLE) 


Fig.  S  —  Error  factors  of  the  optimal  two-point  DDFs  for 
ip  =  2  power  law,  fw  three  different  filter  constraints 
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FREQUENCY  (CYCLES/SAMPLE) 


Fig.  6  —  Error  spectra  after  use  of  the  optimal  two-point 
DDFs,  for  three  different  filtn’  constraints.  The  input  signal 
spectrum  is  a  p  «  2  power  law. 
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FREQUENCY  (CYCLES/SAMPLE) 


Fig.  7  —  Error  factors  of  the  optimal  four-point  DDFs  for 
a  p  =  2  power  law.  for  three  different  filter  constraints 


frequency  (OCLEs/bAMPLE) 


Fig.  8  —  Error  spectra  after  use  of  the  optimal  four-point 
DDFs,  for  three  different  filter  constraints.  The  input  signal 
spectrum  is  a  p  =  2  power  law. 


Notice  that  the  frequency  range  plotted  in  these  figures  extends  to  10  times  the 
Nyquist  frequency.  If  the  power  spectrum  model  of  Eq.  (59)  is  truncated  to  zero  above  the 
Nyquist  frequency,  then,  as  Figs.  3  and  4  show,  the  polynomial  solutions  are  sujwrior  to 
the  "optimal."  Of  course,  the  truncation  produces  a  new  /f(x),  and  a  corresponding  new 
optimal  filter  is  described  by  Eq.  (58).  It  can  be  found  by  a  straightforward  numerical 
procedure  [16],  which  will  not  be  pursued  here.  Instead,  another,  more  robust  approach 
that  achieves  near-optimal  perfOTinance  is  described. 
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4.  SUBOPTIMAL  KERNEL  DESIGN 

Seldom  does  a  simple  analytic  fonn  like  Eq.  (59)  or  (60)  describe  all  power  spectra  of 
interest.  For  imagery,  for  example,  different  power  laws  are  often  fit  to  the  same  spectrum 
in  different  frequency  ranges.  Also,  a  hard  cutoff  at  some  high  frequency  is  generally 
imposed  by  the  optical  system  collecting  the  data.  A  further  complication  is  that  a  spectrum 
may  differ  significantly  in  different  subregions  of  a  single  image.  We  now  discuss  methods 
for  designing  kernels  that  rely  on  more  general  properties  of  power  spectra  than  a  precise 
and  simple  analytic  form. 

4.1  Low-Frequency  Optimal 

The  first  example  is  Symmetric  Low-Frequency  optimal  kernels  of  order  N,  called 
SLF-N.  The  performance  of  these  DDFs  will  be  compared  to  similar  kernels,  LF-N,  found 
for  the  interpolation  approach.  Both  types  of  kernel  are  defined  by  the  requirement  that  the 

corresponding  error  factor  ^^(v)  and  the  maximum  possible  number  of  its  derivatives 

£j'”\v)  be  zero  at  v  =  0.  This  criterion  ensures  a  reasonably  good  performance  for  power 
spectra  with  nearly  singular  low-frequency  behavior,  which  is  common  for  imagery.  The 
LF-N  were  studied  in  Ref.  12;  their  use  is  essentially  a  local  version  of  the  Lagrange 
interpolation  procedure,  which  finds  the  unique  (N-  l)-order  polynomial  connecting  N 
sample  points. 

For  symmetric  DDFs,  Eq.  (53)  implies  that 

4'^{0)  =  0  (66) 

for  all  even  values  of  m.  Therefore,  if  Nc  is  the  number  of  constraints  used  to  preserve  the 
changes  to  be  detected,  then  an  N-point  kernel  has  (N  -  Nc)  degrees  of  freedom,  which  can 
be  used  to  satisfy  Eq.  (66)  for  several  low  odd  values  of  m.  Then  Eq.  (66)  holds  for 
m  =  0,  ....  2{N  ~  Nc).  (In  our  examples,  Nc  always  has  the  value  1.) 

By  contrast,  for  interpolation,  the  N  degrees  of  freedom  can  be  used  to  make  the  error 

in  Eq.  (44)  satisfy  Eq.  (66)  for  all  s  only  for  m  =  0 . N  -  1.  Therefore,  at  low 

frequencies,  the  error  factor  (Eq.  (44))  associated  with  LF-N  varies  as  v^,  while  for  SLF- 
N  with  one  constraint,  it  (Eq.  (43))  varies  as  For  example,  for  N  =  2,  the  LF 

solution  is  just  LIN,  for  which  the  error  factor  varies  as  v^;  for  SLF-2,  derived  below,  it 
varies  as  v^.  Thus,  SLF-N  produces  less  error  at  low  frequencies  than  LF-N.  This 
superior  performance  is  analogous  to  the  situation  found  above  for  the  p  =  2  optimal  DDF, 
although  there  the  DDF  solution  was  even  more  preferred,  in  that  even  the  N  =  2  DDF 
kernel  outperformed  the  optimal  interpolator  for  any  N. 

To  find  the  kernel  for  SLF-N  with  one  constraint  of  the  form  in  Eq.  (56),  we  solve 
Eq.  (66)  which,  from  Eq.  (53),  becomes 

X(n  +  ^)'"ik(n  +  s)  =  0  (67) 

for  m  =  1,  3,  !..,  2N- 3.  For  example,  for  N  =  2,  the  homogeneous  Eq.  (67)  becomes 

|-1  +  ^  /fe(-l  +  ^)  +  2  ® 
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The  inhomogeneous  equation  of  constraint  can  be  taken  to  be  the  dc  condition 

Ik(«+j)=l,  (69) 

ft 

because  the  solutions  Kn  +  for  the  mwe  general  case  of  Eq.  (56)  are  simply  proportional 
to  these.  The  same  is  true  of  the  corresponding  error  factors  Es(v),  which  are  linear  in 
kin  +  s).  The  dc  solution  can  thus  be  used  to  express  the  general  form  of  the  peak-signal  to 
rms  clutter 


SIC='Lk{n+s)t.nlds.  (70) 

ft 

That  is,  if  kin  +  s)  is  the  SLF  solution  with  the  dc  constraint  (Eq.  (69)),  then  Eq.  (7()) 
gives  the  S/C  ratio  for  the  more  general  problem  defined  by  Eq.  (56).  The  denominator  is 
computed  using  Eqs.  (36)  and  (41).  Note  that  the  analogous  property  does  not  generally 
hold  for  the  optimal  kernels  for  a  given  power  spectrum.  Figures  5  and  7  illustrate  this,  in 
that  the  ratio  of  the  error  factors  corresponding  to  different  constraints  clearly  depends  on 
frequency,  whereas  Eq.  (53)  predicts  the  opposite  for  proportional  kernels. 

The  solution  to  Eqs.  (68)  and  (69)  is 

(71) 

This  has  a  simple  interpretation  in  terms  of  conventional  interpolation.  Equation  (71) 
represents  the  kernel  for  linear  interpolation  at  a  shift  of  s/2,  From  Eq.  (52)  then,  the 
filtering  of  signal  2  consists  of  linear  interpolation  at  a  shift  of  -s/2.  In  other  words,  the 
SLF-2  solution  amounts  to  shifting  the  two  discrete  signals  toward  each  other,  using  linear 
interpolation,  by  half  their  separation.  Such  a  simple  relation  does  not,  however,  persist  at 
higher  orders  for  SLF-N. 

The  low-frequency  optimal  dc  kernel  solutions  for  SLF-3  are  straightforward  but 
tedious  to  derive.  They  satisfy  two  homogeneous  equations,  Eq.  (67)  with  m  =  1,  3, 
along  with  Eq.  (69).  The  dc  solution  is 

k(s-l)  = 

k{s)  =-l(s-2)(s  +  2) 

Jk(s+l)  =  ^(s-l)(s-2).  (72) 
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For  N  =  4,  the  dc  solution  is 


^(^-2)=j^(5)(s+  l)(s+2) 

^{5-1)  =  -Y^{5-4)(5+  1)(5  +  2) 

*(^)  =  j^U-4)(s-3)(i+2) 

k{s+l)  =  -^{s-2)is-3){s-4).  (73) 

The  pattern  of  the  monomial  factors  appearing  in  Eqs.  (71)  through  (73)  makes  the 
use  of  an  ansatz  a  feasible  approach  for  higher  iV.  For  example,  when  N=  5,  we  expect 
fourth-order  polynomial  kernels  k(s  -  n)  with  monomial  factors  of  the  form  {s  -  m)  with 
integer  Iml  ^  4,  m  ^  The  coefficients  multiplying  the  various  polynomials  can  be 
determined  from  the  condition  in  Eq.  (69). 

Figure  9  includes  the  error  factors  for  the  above  kernels  along  with  those  for  the 
interpolators  LF-2  (LII^  and  LF-4.  Note  that  the  three-point  DDF  is  better  than  the  four- 
point  inteix)olator  at  all  inband  ftequencies. 


FREQUENCY  (CYCLES/SAMPLE) 

Fig.  9  —  Error  factors  for  two  intopolaiors:  Linear  and  LF-4; 
and  three  DDFs;  SLF-2, 3,  and4 

Figures  10  and  1 1  show  the  SLF-2,  3,  and  4  error  spectra  for  p  =  2,  4  power  laws, 
respectively.  Note  that  no  meaningful  comparison  can  be  made  between  Figs.  10  and  1 1 
because  no  low-frequency  cutoffs  have  been  defuied  that  could  m^e  the  input  spectra  have 
equal  areas  and,  hence,  equal  (and  finite)  unfiltered  variances.  Forp  =  2,  the  rms  difference 
ds  is  approximately  3.3  times  smaller  for  SLF-3  than  for  SLF-2;  and  7.1  times  smaller  for 
SLF-4  than  for  SLF-2.  For  p  =  4,  they  are  3.8  and  8.7.  (For  p  =  3,  the  corresponding 
ratios  are  3.5  and  7.7,  respectively.)  The  more  singular  spectra  benefit  relatively  more  as  N 
increases  for  the  SLF  filters,  which  are  designed  to  work  well  at  low  frequencies. 
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Fig.  10  —  Error  spectra  for  SLF-2, 3,  and  4, 
and  a  p  =  2  power  law 


FREQUENC  i'  (CYCLES/SAMPLE) 


Fig.  11  —  Error  spectra  for  SLF-2, 3,  and  4, 
and  ap  =  4  power  law 


Notice  that  the  residual  error  spectrum  plotted  in  Fig.  10  varies  as  near  dc,  and 

as  in  Fig.  11.  Even  foriV  =  2,  this  may  be  considered  low-frequency  overkill.  The 
bulk  of  the  error  energy  even  for  the  p  =  4  spectrum  resides  at  high  frequencies,  which 
have  been  neglected  in  the  SLF  definition.  The  next  subsection  expands  the  filter  design 
procedures  to  remedy  this  defect. 
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4.2  Hybrid  Designs 

The  error  factor  in  Eq.  (53)  and  derivatives  of  it  may  be  constrained  to  the  value  zero 
at  other  frequencies  besides  dc.  For  many  of  the  DDFs  and  interpolators  already  discussed, 
the  peak  values  of  the  error  spectra  occur  at  a  high  frequency,  in  particular,  above  the 
Nyquist.  We  now  restrict  the  remaining  discussion  to  oversampled  systems,  such  as  are 
assumed  in  Figs.  9  through  11.  For  them,  the  maximum  error  occurs  at  the  Nyquist 
frequency,  and  a  natural  next  step  in  filter  design  would  be  to  sacrifice  some  low-frequency 
performance  to  reduce  the  residual  error  at  higher  frequencies. 

The  elimination  of  errors  at  nonzero  frequencies  is  possible  with  interpolation  as  well 
as  with  DDFs  and  has  been  studied  [12].  Notice,  however,  that  Es(y)  for  interpolation  (Eq. 
(44))  is  complex  at  most  frequencies,  and  requiring  it  to  be  zero  usually  expends  two 
degrees  of  design  freedom.  An  exception  occurs  at  dc,  at  which  frequency  requiring  Eq. 
(44)  to  be  zero  at  all  shifts  s  constitutes  but  one  real  constraint.  On  the  other  hand,  because 

at  V  =  1/2  the  transfer  function  //j*^  of  Eq.  (42a)  is  real,  £,(v)  in  Eq.  (44)  (i.e.,  for 
interpolation)  can  never  be  zero  at  all  shifts.  This  limitation  occurs  because  a  sine  wave  at 
the  Nyquist  frequency  is  zero  at  all  integer  sample  values,  and  is,  therefore,  invisible  to  the 
sampling  process.  Such  a  sinusoid  cannot,  therefore,  be  reproduced  by  any  interpolator. 

On  the  other  hand,  for  symmetric  DDFs,  Eg(\)  may  be  taken  as  real  (see  Eq.  (53)). 
Therefore,  each  frequency  null  requires  the  expendimre  of  only  one  degree  of  freedom. 
Furthermore,  because  of  (52),  the  transfer  functions  in  Eq.  (42)  for  the  two  filters  are 
complex  conjugates,  and  the  complex  error  factor  in  Eq.  (43)  may  be  written  as 

£,(v)  -4  2|w{'>(v)|sin(()><'>(v)-^) ,  (74) 

in  which  <|>(^)  is  the  phase  of  A  unimodular  factor  has  again  been  (Knitted. 


Now  note  from  Eq.  (45)  that  if  in  Eq.  (74)  is  derived  from  the  sine  kernel  at  the 

shift  s/2,  then  <|>  has  the  value  oar/2  inband,  making  E/y)  zero  inband  for  all  shifts.  This 
simply  confirms  that  shifting  a  pair  of  discrete  oversampled  signals  toward  each  other  by 
half  their  separation  gives  the  ideal  result  if  the  ideal  interpolator  can  be  used.  However, 
Eq.  (74)  shows  that  zero  inband  error  for  DDFs  does  not  require  use  of  the  ideal 
interpolator,  but  only  one  with  the  ideal  phase.  By  contrast,  for  the  interpolation  problem, 
zero  inband  error  in  Eq.  (44)  requires  specific  use  of  the  sine  interpolator. 

According  to  Eq.  (74),  E/y)  can  be  made  zero  either  by  making  the  magnitude  of  the 
transfer  function  zero,  or  by  adjusting  its  phase  appropriately.  Furthermore,  this  can  be 
accomplished  with  DDFs  even  at  the  Nyquist  frequency,  an  impossibility  for  interpolators. 

As  an  example,  we  consider  three-point  DDFs  for  which  the  error  is  zero  at  angular 
frequencies  coi  and  (02.  The  form  in  Eqs.  (53)  and  (56)  can  be  written  as 

X\  x/ 

v2 

A.i  Aq  Aj 

.  1  t-i  . 


^(-1+  s) 

o' 

k(s) 

= 

0 

k(\+s) 

.  1  . 

(75) 
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where 

4  =  sin[(0,(n+i)]  (i  =  l.:  /i  =  -l,0.  1).  (76) 

According  to  Cramer's  Rule,  the  solution  kin  +  s)  can  be  derived  by  substituting  the  right- 
hand  vector  in  Eq.  (75)  for  the  nth  column  of  the  left-hand  matrix,  finding  the  determinant 
of  this  modified  matrix,  and  dividing  by  the  determinant  o'"  the  original  matrix.  However, 
the  zeroes  in  the  vector  make  each  such  modified  matrix  independent  of  all  Therefore, 
any  solutions  [kin  +  5)}  for  one  set  (/„)  are  proportional  to  those  for  any  other,  for 
example  for  f  /„  =  1 } .  We  may,  therefore,  solve  Eq.  (75)  for  /„  =  1  and  rely  on  Eq.  (70)  for 
the  physically  important  quantity,  peak-signal/clutter.  Note  that  because  the  SLF-3  filter 
may  be  viewed  as  the  limiting  solution  to  Eq.  (75)  as  the  frequencies  Vj,  V2  0,  the 
prescription  of  Eq.  (70)  also  applies  to  it  and,  more  generally,  to  SLF-N. 

The  solution  to  Eq.  (75)  with  =  1  is 

k{s-l)  =  mii  Idet 
kis)  =  -m32  I det 

ik(5+ 1)  =  W33 /^er  (77) 

with  m,)  the  cofactors  of  the  matrix  in  Eq.  (75)  ariu  aet  its  determinant 

Figure  12  plots  the  error  factor  for  the  above  kernel  with  Vi  =  1/4,  V2  =  1/2.  We  will 
generally  use  the  notation  5//(Vi,  V2, ...)  to  mean  a  symmetric  A/-point  DDF  with  zeroes  in 
its  error  factor  at  frequencies  Vi,  V2, ... .  If  the  number  of  arguments  of  SN  is  less  than 
A  -  1  —  the  single  dc  constraint  of  Eq.  (69)  is  always  assumed  —  then  the  remaining 
degrees  of  freedom  are  understood  to  have  b^n  used  to  set  odd  low-order  derivatives  of 
^s(v)  equal  to  zero  at  v  =  0.  (Even  ones  are  automatically  zero.)  Thus,  53(1/4),  whose 
performance  is  also  depicted  in  Fig.  12,  satisfies  £,(1/4)  =  0  as  well  as  £,(0)  =  0.  Such 
hybrid  filters  generally  have  low-frequency  properties  similar  to  SLFs,  along  with  null 
errors  at  selected  higher  frequencies.  As  with  SLFs,  for  hybrid  filters  the  dc  constraint  may 
be  assumed  and  Eq.  (70)  used  to  find  S/C  ratios. 

Figure  12  also  shows  the  performance  of  an  A  =  4  hybrid  interpolator  studied  in  Ref. 
12.  In  this  case,  two  of  the  inte^olator's  four  degrees  of  freedom  are  used  to  create  the  null 
at  V  =  1/4;  the  other  two  make  £,(0)  =  £i(0)  =  0.  Therefore,  at  low  v,  l£,l  varies  as  v^  for 
the  interpolator.  For  53(1/4,  1/2),  l£,l  varies  as  v,  performing  worse  than  the  interpolator  at 
low  frequencies,  but  better  at  all  frequencies  beyond  v  «  0.13.  The  error  factor  for  53(1/4) 
begins  with  a  dependence  and  remains  superior  to  the  interpolator  at  all  frequencies,  but 
becomes  worse  than  the  other  DDF  above  v  «  0.39. 

Figure  13  plots  the  corresponding  error  spectra  for  the  p  =  2  power  spectrum.  Both 
DDFs  give  smaller  residual  rms  differences  than  the  interpolator  (53(1/4, 1/2)  by  the  factor 
3.4  and  53(1/4)  by  the  factor  4.5). 
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FREQUENCY  (CYCLES/SAMPLE) 


Fig.  12  —  Errors  factors  for  three  hybrid  methods.  Besides  at  the 
frequency  v  =  0,  other  nulls  occur  at  1/4  for  the  interpolator  and  for 
both  DDFs,  S3(l/4, 1/2)  and  S2(l/4);  also  at  v  =  1/2  for  S3(l/4. 1/2). 
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FREOUENC')  (CYCLES/SAMPLE) 

Fig.  13  —  Error  spectra  for  three  hybrid  methods  and  ap-2  spectrum. 
The  interpolator  is  of  order  4;  both  DDFs  are  of  order  3. 


The  relative  gain  in  performance  is  even  more  impressive  if  the  order  of  the  DDF  is  as 
high  as  that  of  the  interpolator.  Figure  14  compares  the  error  factors  for  S4(l/4,  1/2)  with 
the  two  DDFs  of  Fig.  12.  The  extra  degree  of  freedom  in  the  four-point  method  allows  it 
to  combine  the  low-frequency  behavior  of  53(1/4)  with  the  high-frequency  behavior  of 
53(1/4, 1/2).  The  corresponding  error  spectra  are  shown  in  Fig.  15.  Figure  16  summarizes 
the  /7  =  2  results  and  includes  results  for  53(0.15, 0.4)  and  54(0.1,  0.3,  0.45),  which  have 
been  found  to  be  good  general  purpose  DDFs  for  imaging  applications.  Similar  gains  are 
found  for  p  =  3  and  p  =  4  power-law  spectra. 
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FREOUEUS  :  (C'l'CLES/'SAMPLEl 

Fig.  14  —  Comparison  of  hybrid  DDFs 


FREQUENCY  (CYCLES/SAMPLE) 


Fig.  15  —  Error  spectra  for  hybrid  DDFs  and  ap  =  2  input  spectrum 
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Fig.  16  —  Performance  of  interpolators  and  DDFs  relative  to  Linear, 
for  a  /7  =  2  input  spectrum  truncated  at  the  Nyquist  frequency.  Numb^ 
in  parentheses  indicate  locations  of  null  frequencies. 


4.3  Effects  on  Gaussian  Targets 

Figures  17  through  20  illustrate  the  effects  of  two  four-point  DDFs  on  Gaussian 
targets.  Figures  17  and  18  use  the  p  =  2  optimal  DDF,  with  Gaussian  target  constraints  in 
Eqs.  (56)  and  (57)  for  o  =  0.5  and  1.0,  respectively.  Figures  19  and  20  show  similar 
results  for  the  DDF  54(1/4,1/2).  The  targets  are  assumed  to  have  moved  a  distance  of  at 
least  a  few  samples  between  the  collection  times  of  the  discrete  frames  of  data,  so  that  they 
are  not  destroyed  by  the  frame  subtraction  operation.  The  solid  lines  represent  the 
Gaussians,  which  are  sampled  at  points  on  the  abscissa  separated  by  a  distance  of  one.  If, 
for  example,  a  Gaussian  is  sampled  at  its  peak,  located  at  position  5  in  the  plots,  then  the 
filtered  target  has  a  peak  value  identical  to  the  unfiltered,  t^ause  the  two  curves  intersect 
there,  a  consequence  of  imposing  the  condition  of  (56)  on  the  filters.  If,  on  the  other 
hand,  a  Gaussian  signal  happens  to  be  sampled  off  its  peak,  then  the  corresponding  value 
produced  by  the  filter  is  usually  higher  than  the  input.  That  is,  the  filtered  curves  are 
usually  above  the  unfiltered  at  any  given  sampling  position. 


27 


AMPLITUDE  AMPLITUDE 


A.SCHAUM 


Fig,  17  —  Effect  on  a  o  =  0,5  Gaussian  target  of  the  optintal 
four-point  DDF,  designed  with  a  a  »  0.S  Gaussian  constraint, 
for  a  p  a  2  power-law  input  spectrum  and  the  shift  s  -  0.25 
samples 


Fig.  18  —  Effect  on  a  o  a  1.0  Gaussian  target  of  the  optimal 
four-point  DDF,  designed  with  a  a  a  1.0  Gaussian  constraint, 
for  a  p  a  2  power-law  input  spectrum  and  the  shift  s  =  0.25 
samples 
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Fig.  19  —  Effect  of  54(1/4, 1/2)  DDF  with  a  =  O.S  Gaussian 
constraint  on  a  o  =  03  Gaussian  target 


POSITION  (SAMPLES) 

Fig.  20  —  Effect  of  54(1/4, 1/2)  DDF  with  o  =  1.0  Gaussian 
constraint  on  a  a  =  1.0  Gaussian  target 


Thus,  the  signal-to-clutter  formula  (Eq.  (70))  is  usually  a  slight  underestimate  of  the 
true  gain.  For  example,  if  in  Fig.  17  the  sampling  points  arc  at  the  half  integers,  then  the 
peak  input  signal  is  exp(-(l/2)2/(2o2))  =  0.607,  whereas  the  peak  filtered  value 
is  =  0.75.  Therefore,  Eq.  (70)  actually  underestimates  the  peak-signaVclutter  by  the  factor 
1.23.  Note  also  that  a  second,  filtered  version  of  a  moving  Gaussian  target  appears  in  the 
difference  frame;  for  symmetric  DDEs,  it  corresponds  to  the  negative  mirror  image  of  the 
appropriate  dotted  curve  in  Figs.  17  through  20. 
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More  generally,  the  perfcamance  of  the  single-frame  filters  can  be  characterized  by  the 
usual  frequency  descriptor  of  digital  signal  processing,  the  magnitude  of  the  transfer 
function  Figure  21  plots  these  for  three  three-point  DDFs  with  the  dc  constraint. 
These  show  the  typical  low-pass  characteristic  of  all  the  filters  discussed  here. 


FREQUENCY  (CYCLES/SAMPlE) 


Fig.  21  —  Transfer  funcdcms  of  three  three-point 
DDFs  with  dc  constraint 

Figure  22  shows,  for  each  of  the  three  constraints,  the  transfer  functions  for  an 
effective  inband  DDF,  the  hybrid  54(1/4,1/2).  Note  that  four  of  the  six  filters  described  in 
Figs.  21  and  22  have  difference-error  nulls  at  v  =  1/2,  five  have  them  at  v  =  1/4,  and  all 
have  them  at  v  =  0.  However,  these  nulls  are  achieved  differently,  according  to  Eq.  (74). 
The  null  error  factors  at  v  =  1/2  are  always  zero  because  I//W|  =  0.  By  contrast,  the  plots 
show  that  is  never  zero  at  the  two  lower  frequencies,  and  so  these  nulls  are  achieved 
by  having  =  (Os/2. 


Note  also  that  the  three  transfer  functions  in  Fig.  22  are  proportional,  a  reflection  of 
the  previously  described  proportionality  of  the  corresponding  kernels  {k(/z  +  5)}  for  hybrid 
DDFs.  This  should  be  contrasted  with  the  optimal  four-point  filter  for  the  p  =  2  spectrum  in 
Fig.  23,  where  the  dependence  on  the  inhomogeneous  target  constraint  of  Eq.  (56)  is  more 
complicated.  Rnally,  in  each  of  these  two  figures,  the  two  dc  amplifications  occur  because 
of  the  combination  of  the  low-pass  filter  characteristics  of  //(*>  and  the  constraint  preserving 
peak  value  of  a  Gaussian  target.  These  require  the  filter  to  broaden  the  target  without 
depressing  its  central  value.  The  dotted  curves  in  Figs.  17  through  20  all  have  larger 
means  than  the  Gaussians  that  generated  them. 

5.  RELATED  ISSUES 

Reference  16  discusses  several  corollaries  of  the  mean  squared-error  formulation 
developed  here.  These  include  generalizations  to  higher  dimensiems,  in  which  the  variables 
X,  n,  and  the  shift  5  become  vectors.  Also,  an  extension  of  the  analysis  to  include  more 
than  two  signals  is  achieved  which,  for  image  processing  applications,  transforms  the 
problem  into  one  of  multiframe  tracking.  Geno^izations  in  Ref  16  are  also  developed  that 
allow  the  evaluation  and  optimization  of  methods  of  prediction,  restoration,  and  finite- 
length  matched  filtering. 
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Fig.  22  —  Transfer  functions  of  S4(l/4, 1/2)  for 
three  diffetent  target  constraints 


FREQUENCY  (CYCLES/SAMPLE) 


Fig.  23  —  Magnitude  of  transfer  functions  of  the  four-point 
optimal  (for  a  p  =  2  input  power-law  spectrum)  DDFs 
corresponding  to  three  different  target  constraints 
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Here  we  include  a  discussion  of  additive  temporal  noise.  This  report  is  concerned 
mostly  with  the  error  dg,  defined  in  Eq.  (25).  In  image  processing  applications,  dg  is  called 
residual  clutter.  Initial  clutter  is  just  the  variance  of  the  unprocess^  image.  According  to 

Eqs.  (41)  and  (36),  if  the  second  filter  k2  vanishes  and  K^^\n)  is  5q,  i.e.,  if  nothing  is 

done  to  the  first  image,  then  1  and  dj  is  the  variance  of  the  image  (by  Plancherel's 
formula  [19]).  Therefore,  according  to  Eq.  (36),  the  squared  modulus  of  the  error  factor  is 
the  ratio  of  dte  power  in  the  filtered  difference  spectrum  to  that  in  the  original  image. 

In  images  for  which  stochastic  (temporal)  noise  appearing  in  both  images  is 
comparable  to  the  spatial  clutter,  or  for  processed  difference  images  in  which  the  residual 
clutter  has  been  reduced  to  the  noise  level,  Eq.  (36)  is  inaccurate.  Stochastic  noise 

contributions  M  must  be  added  to  the  digital  signals,  so  that  the  replacements 

f{n  +  r)  ->/(n  +  s)  +  N2(«)  (78) 


are  required. 

Note  that  if  the  second  image  is  unaltered,  in  which  case  dg  is  a  measure  of 
interpolation  accuracy,  then  the  methodology  of  Sections  2  and  3  applies  to  the  restoration 

of  a  digital  image  with  additive  noise  N.  In  this  case,  the  problem  becomes  one  of  finite 
Wiener  filtering,  and  the  optimization  equation  describes  the  minimum  mean-squared  error 
finite  digital  filter.  Note  that  this  is  usually  not  simply  a  truncated  version  of  the  standard 
Wiener  solution,  which  involves  an  infinite-order  filter  kin  +  s). 

Here  we  assume  that  the  noise  terms  in  Eq.  (78)  are  independent  both  of  each  other 
and  of  the  image/,  and  that  they  are  zero-mean  stationary  ergc^c  random  processes  with 

the  following  common  autocorrelation  function  A.: 

Ain)  =  Zl?^i{m)f<i{n-r4 ,  (79) 

where  £  denotes  an  ensemble  average  for  i  =  1  and  2.  Then  Eq.  (54)  generalizes  to 

<<?  =  2  S  ([«(»)+ .4  (n)]t(n  +  m  +  s)k{m  +  s)- R{n  +  s) k[n  +m  +  s) s)) ,  (80) 
m,n 

which  includes  both  clutter  and  noise.  The  optimization  (Eq.  (55a))  is  replaced  by 


'Z[R{m)  +  A{m^k(n -m  + s)-'ZR{m  + s)k{m -n  +  (81) 

m  m  i  *  1'*+'^) 
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As  an  example  of  differences  between  optimal  solutions  with  and  without  noise,  we 
consider  N  =1  solutions  for  the  Lorentzian  clutter  spectrum,  with  the  general  target 

constraint  of  Eq.  (56)  (with  to  =  1).  The  temporal  noise  is  assumed  to  be  white,  i.e.,  A.{m) 

oc  5o ,  and  r  is  the  noise  to  clutter  ratio  r  =  #4.(0)//?(0).  The  optimal  DDF  is  given  by 


^(-1) 


_ p-  t\  (p^-1-  r) _ 

p2-j _  1_ _  2ti  (p -  pl-^  -  rf(l-  p^  +  r) 


/t(0)=l-ri*(-l). 
with  p  given  by  Eq.  (61). 

If  the  noise  were  absent,  then  r  =  0  and  the  solution  as  5  ->  0  becomes 


(82) 


*(-!)  =  0 


*(0)=1.  (83) 

This  corresponds  to  complete  removal  of  clutter  by  the  DDF  by  reproducing  the  sample 
values  of  both  digital  signals.  If  the  solution  in  E!q.  (82)  with  r  =  0  is  used  to  reduce  a 
severe  clutter  problem  that  might  exist  at  large  shifts — that  is,  if  noise  is  ignored — then  at 
small  shifts  it  can  be  noticeably  subopdmal,  as  measured  by  the  ratio  of  signal  to  clutter- 
plus-noise.  In  particular,  if  the  noise  term  is  kept  in  Eq.  (82),  the  solution  at  5  =  0  is 


k  (-1)  = _ ISx _ 

l-p2  +  r(l-hrf) 

k(0)  =  l -tiki-l) . 


(84) 


Noise  cannot  be  totally  removed  at  any  shift,  and  this  solution  represents  a  balance  in 
which  the  filter  keeps  some  clutter  and  some  noise.  If  noise  happens  to  dominate  clutter, 
i.e.,  r  is  large,  then  Eq.  (84)  approaches 


«-l) - '-hr 

l+'i 

m — 1— ,  (85) 

l+(f 

which  is  the  two-point  matched  filter  in  white  noise.  This  DDF  produces  a  smaller  (mean 
squared)  signal-to-noise  ratio  than  does  the  DDF  of  Eq.  (83)  by  the  factor  (1  +  Such  is 
the  extent  to  which  the  DDF  designed  purely  for  clutter  suppression  is  suboptimal,  when  it 
is  applied  in  a  noise-dominant  environment. 

6.  SUMMARY  AND  CONCLUSIONS 

This  report  generalizes  the  notion  of  interpolation  an  J  digital  subtraction  as  a  means  of 
detecting  changes  in  a  pair  of  discrete  signals.  By  acting  on  both  signals,  Dual  Difference 
Filters  afford  more  flexibility  and  power  than  do  interpolators  of  the  same  order.  This  is 
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proven  by  a  fundamental  theorem  (Eq.  (36))  of  mean  squared  error  for  the  general 
resampling  problem,  which  includes  inteipolaticm  as  a  special  case. 

For  any  Hxed  shift  between  sampling  grids,  the  strong  form  of  the  theorem  associates 
an  error  spectrum  with  any  DDF  that  is  applied  to  an  oversampled  signal.  The  error 
spectrum  is  the  product  of  the  power  spectrum  of  the  continuous  underlying  signal  that 
generates  the  discrete  samples,  and  an  attenuation  factor.  The  weak  form  of  the  theorem 
says  the  same  thing  for  undersampled  signals,  but  in  an  average  sense.  Formulae  are 
derived  for  the  attenuation  factors  that  can  be  used  for  filters  that  are  of  finite  extent  in  either 
frequency-space  (Eq.  (37))  or  real-space  (Eq.  (41)).  A  corollary  of  the  weak  theorem 
allows  application  of  the  error  fonnulae  to  stochastic  problems. 

A  fundamental  optimization  equation  (Eq.  (58))  is  derived.  If  the  precise  form  of  the 
signal  spectrum  is  known,  then  the  solution  of  tiie  equation  produces  the  maximum 
possible  peak-signal  to  rms  clutter  for  any  DDF  of  given  order.  As  an  example,  the  case  of 
a  power-law  spectrum  and  a  single  linear  constraint  on  the  filters  is  solved. 

Flexible  design  methods  for  DDFs  are  developed  that  produce  superior  performance 
for  a  wide  range  of  signal  spectra.  For  example,  for  three  bandlimited  power-law  spectra 
(p  =  2, 3,  and  4),  a  single  four-point  DDF  is  found  that  enhances  residual  signal/clutter  by 
approximately  30  dB.  compared  to  the  best  popular  interpolates  of  the  same  eder. 
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